#include "yocto/math/dat/spline.hpp"
#include "yocto/math/io/data-set.hpp"
#include "yocto/ios/icstream.hpp"
#include "yocto/ios/ocstream.hpp"

#include "yocto/exception.hpp"
#include "yocto/code/hsort.hpp"
#include "yocto/string/conv.hpp"
#include "yocto/math/fit/lsf.hpp"
#include "yocto/math/dat/linear.hpp"
#include "yocto/string/vfs-utils.hpp"

using namespace yocto;
using namespace math;

class SplineFit
{
public:
    explicit SplineFit() : 
    Xa(),
    BoundLo(false,0),
    BoundHi(true,0)
    {
    }
    
    virtual ~SplineFit() throw()
    {
    }
    
    double Eval( double x, const array<double> &Ya )
    {
        assert( Ya.size() == Xa.size() );
        spline<double> s(Xa,Ya,BoundLo,BoundHi);
        return s(x);
    }
    
    bool Callback(const fit::sample<double>              &s, 
                  const fit::lsf<double>::field          &,
                  const array<double>                    &Ya )
    {
        std::cerr << "D=" << s.D << std::endl;
        std::cerr << "Ya=" << Ya << std::endl;
        return true;
    }
    
    vector<double> Xa;
    const spline<double>::boundary BoundLo;
    const spline<double>::boundary BoundHi;
    
private:
    YOCTO_DISABLE_COPY_AND_ASSIGN(SplineFit);
    
};

int main( int argc, char *argv[] )
{
    
    try 
    {
        
        if( argc <= 2 )
            throw exception("need input.dat xcontrol.dat");
        
        const string input_file(argv[1]);
        const string xctrl_file(argv[2]);
        
        //-- load sample
        vector<double> X;
        vector<double> Y;
        {
            ios::icstream fp( input_file );
            data_set<double> ds;
            ds.use(1,X);
            ds.use(2,Y);
            ds.load( fp );
        }
        const size_t        N = X.size();
        if( N < 2 )
            throw exception("Not enough data points");
        vector<double>      Z(N,0.0);
        co_hsort(X, Y);
        
        //-- create fit sample
        fit::sample<double> S( X, Y, Z );
        
        //-- load control points
        SplineFit Sf;
        vector<double> &Xa =Sf.Xa;
        {
            ios::icstream fp( xctrl_file );
            data_set<double> ds;
            ds.use(1,Xa);
            ds.load(fp);
        }
        
        Xa.push_back( X[1] );
        Xa.push_back( X[N] );
        
        hsort(Sf.Xa);
        
        //-- check
        const size_t   Na = Xa.size();
        for( size_t i=1; i<=Na; ++i )
        {
            if( Xa[i] < X[1] || Xa[i] > X[N] )
                throw exception("invalid Xa=%g\n", Xa[i]);
        }
        
        //-- initialize from data points
        vector<double> Ya(Na,0.0);
        for( size_t i=1; i <= Na; ++i )
        {
            Ya[i] = linear(Xa[i], X, Y);
        }
        vector<bool>   used(Na,true);
        vector<double> aerr(Na,0.0);
        used[1]  = false;
        used[Na] = false;
        
        //-- fit
        fit::lsf<double>::field     F( &Sf, & SplineFit::Eval );
        fit::lsf<double>::callback  Cb( &Sf, &SplineFit::Callback);
        fit::lsf<double>        LeastSquare;
        std::cerr << "Fitting..." << std::endl;
        std::cerr << "Ya=" << Ya << std::endl;
        LeastSquare.ftol = 1e-4;
        LeastSquare( S, F, Ya, used, aerr, &Cb );
        std::cerr << "Ya=" << Ya << std::endl;
        
        {
            string  outfile = input_file;
            _vfs::change_extension(outfile, "spline.dat" );
            ios::ocstream fp(outfile,false);
            for( size_t i=1; i <= N; ++i )
            {
                fp("%g %g %g\n", X[i], Z[i], Y[i] );
            }
        }
        return 0;
    } 
    catch( const exception &e )
    {
        std::cerr << e.what() << std::endl;
        std::cerr << e.when() << std::endl;
    }
    catch (...)
    {
        
    }
    return 1;
}
